%% Stochastic Simulations in Nonlinear Models
% by Jaromir Benes
%
% Set up and run simulations of unanticipated shocks. Non-linear
% simulations in IRIS are equivalent to a perfect-foresight solution. This
% does not mean unanticipated stochastic shocks cannot be simulated: it
% only means that simulations of such shocks must be split into a number of
% overlapping sub-simulations (called segments in IRIS), depending on the
% occurence of unanticipated shocks.

%% Clear Workspace

clear;
close all;
home;
irisrequired 20140401;
%#ok<*NOPTS>

%% Load Model Object

load read_model m;

%% Parameterize Std Deviations of Shocks
%
% Assign some values to the std deviations of shocks. These values are used
% when shocks are randomly drawn and simulated. There is no need to resolve
% the model after new std deviations have been changed since they do not
% affect the first-order approximate solution matrices.
%
% Keep the std deviations rather low, otherwise the nonlinear simulations
% may not converge for very low levels of monetary policy credibility.

m.std_er = 0.2;
m.std_ey = 0.4;
m.std_epi = 0.4;
m.std_et = 0;

%% Create Steady-State Database with Random Shocks
%
% Call the function `sstatedb` to create a database with each model being
% variable assigned its steady-state value. The range on which these time
% series are created on a range -2:20, and not just 1:20 as the second
% input argument says. This is because `sstatedb` automatically looks up
% the maximum lag in the model and adjusts the range accordingly to include
% all necessary initial conditions so that a simulation can start at time
% t=1.
%
% With the option `'randomShocks=' true` <?randomShocks?>, the model shocks
% are randomly drawn on the range 1:20 from a normal distribution with mean
% zero and std deviations given by the model object.

rng(0);
disp('Currently assigned std deviations of shocks')
disp(get(m,'std'));

d = sstatedb(m,1:20, ...
    'randomShocks=',true) %?randomShocks?

d.y

%% Linearized Simulation
%
% First, simulate the random shocks using first-order approximate solution
% (linearized model). The shocks are simulated as unanticipated. Note that
% credibility plays no role in this simulation, because the level of
% credibility has no first-order effects in the model equations.

d1 = d;
s1 = simulate(m,d,1:20,'dbOverlay=',true,'anticipate=',false);

%% Noninear Simulations
%
% Next, simulate the random shocks in an exact nonlinear mode. Because
% the shocks are unanticipated, the whole simulation range, 1:20, must be
% segmented by the occurence of unanticipated shocks. As a result, a total
% of 20 overlapping sub-simulations will be run. It works as follows:
%
% * Because there is an unanticipated shock in every period, the simulation
% range will be divided into 20 segments, each one period long.
%
% * The simulation of the 1st segment will take initial conditions from
% the input database, `d2` or `d3`, and the shock in period t=1. All other
% shocks (t=2, t=3, etc.) will be disregarded. The nonlinear simulation of
% the 1st segment will be performed 15 periods long (this is controlled
% by the option `'nonlinearise='`). However, only the first simulated
% period, i.e. t=1, will be stored and included in the output database.
%
% * The simulation of the 2nd segment will take initial conditions from
% the the simulation of the first segment, and the shock in period t=2 from
% the input database. Again, the 2nd segment will be simulated for
% another 15 periods (i.e. t=2 through t=16), but only the first simulated
% period, i.e. t=2, will be included in the output database.
%
% * The simulation of the third and all other segments will follow in the
% same way.
%
% * The output database will include period t=1 from the simulation of the
% 1st segment, period t=2 from the simulation of the 2nd segment, ...
% period t=20 from the simulation of the 20th segment.
%
% * Within each segment, the simulation takes several iterations. These
% iterations are, by default, reported on the screen: A one-line report is
% printed for every 100 iterations. The first piece of information at each
% line is the segment in the following format:
%
%    XX:YY[ZZ]#NN
%
% where `XX` is the first period of the respective segment, `YY` is the
% last period stored for the output database from this segment (i.e.
% periods `XX` : `YY` will be included in the output database), `ZZ` is the
% actual last period simulated in a nonlinear mode, and `NN` is the total
% number of simulated periods. If `ZZ` > `YY`, then additional simulation
% periods are added, if `ZZ` < `YY`, then a total of `YY` periods are
% simulated, but only `ZZ` periods are in a nonlinear mode, and all
% periods beyond that, i.e. `ZZ`+1 : `YY` are simulated using first-order
% solution.
%
% The above segmented simulation can be depicted as follows:
%
% <latex>
% \medskip
% \renewcommand{\arraystretch}{1.1}\renewcommand{\tabcolsep}{0.5em}
% \newcommand{\mycol}[1]{\makebox[1em]{#1}}
% \newcommand{\mysgm}[1]{\ttfamily\small#1}
% \begin{tabular}{|r@{ }l|c|c|c|c|c|c|c|c|c|c|c|c|c|c|c|c|c|c|}
% \hline
% & \hfill $t=$ & \mycol{1} & \mycol{2} & \mycol{3} & \mycol{4} & \mycol{5} & \mycol{6} & \mycol{7} & \mycol{8} & \mycol{9} & \mycol{10} & \mycol{11} & \mycol{12} & \mycol{13} & \mycol{14} & \mycol{15} & \mycol{16} & \mycol{17} & \mycol{18} \\ \hline
% 1st & segment \mysgm{1:1[15]} & X & x & x & x & x & x & x & x & x & x & x & x & x & x & x &   &   &   \\ \hline
% 2nd & segment \mysgm{2:2[16]}&   & X & x & x & x & x & x & x & x & x & x & x & x & x & x & x &   &   \\ \hline
% 3rd & segment \mysgm{3:3[17]}&   &   & X & x & x & x & x & x & x & x & x & x & x & x & x & x & x &   \\ \hline
% \multicolumn{20}{|c|}{$\vdots$} \\
% \multicolumn{20}{|c|}{$\vdots$} \\ \hline
% \end{tabular}
% \medskip\par
% where an X\textit{s} denote periods that will be included in the output
% database, whereas x\textit{s} denote periods simulated in that particular
% segment, but not included in the output database.
% </latex>


% ...
%
% Simulate random shocks with full initial credibility, `c(0) = 1`.

d2 = d;
d2.c(0) = 1;
s2 = simulate(m,d2,1:20,'dbOverlay=',true,'anticipate=',false, ...
    'nonlinear=',15,'maxiter=',5000);

% ...
%
% Simulate random shocks wih low initial credibility, `c(0) = 0.3`.

d3 = d;
d3.c(0) = 0.3;
s3 = simulate(m,d3,1:20,'dbOverlay=',true,'anticipate=',false, ...
    'nonlinear=',20,'maxiter=',5000);

%% Plot Simulation Results
%
% Plot the simulated paths for all three simulations in one graph to
% compare the effect of nonlinearities:
%
% * The level of credibility does not move throughout the linear simulation
% (and even though it had moved it would have had no effect on the rest of
% the model anyway). This is very different in the other two simulations
% where credibility is affected by inflation outcomes.
%
% * Because of reductions in credibility in both of the nonlinear
% simulations, it costs more output to bring inflation back to the target.
% This is because the Phillips curve becomes more backward-looking with
% low levels of credibility
%
% * The other factor making inflation stabilisation costlier is the
% asymmetry of the Phillips curve in the output gap itself. Higher levels
% of inflation require a deeper slowdown in real economic activity (and the
% same positive shocks to the output gap will cause inflation to shoot up
% more).

dbplot(s1 & s2 & s3,1:20, ...
    {'"Inflation" pi', ...
    '"Credibility" c', ...
    '"Output gap" y', ...
    '"Policy rate" r'}, ...
    'zeroline=',true,'tight=',true,'marker=','.');

grfun.bottomlegend('Linearized Simulation', ...
    'Nonlin Simulation with High Initial Credibility', ...
    'Nonlin Simulation with Low Initial Credibility');

%% Inverted Simulation
%
% Use the above simulated database, `s3`, as an input database, and swap
% the endogeneity and exogeneity of variables and shocks on the entire
% simulation range. In other words, exogenize (fix) the simulated paths for
% the endogenous variables, `y`, `pi`, `r` <?exogenize?>, while
% endogenizing the shocks `ey`, `epi`, `er` <?endogenize?>. The function
% `simulate` uses only the information on initial condition and the
% exogenized data points from the input database. Remember to enter the
% information on endogenized and exogenized data points to the function
% `simulate` through the option `'plan='` <?plan?>.
%
% The simulated paths for the endogenized shocks, `ey`, `epi`, `er`, in the
% output database `s4` are, of course, identical to the shocks randomly
% generated earlier with numerical precision corresponding to the tolerance
% level in nonlinear simulations (which is 1e-5 by default, but can be
% changed through the option `'tolerance='`);

p = plan(m,1:20);
p = exogenise(p,{'y','pi','r'},1:20); %?exogenize?
p = endogenise(p,{'ey','epi','er'},1:20); %?endogenize?

s4 = simulate(m,s3,1:20,'dbOverlay=',true,'anticipate=',false, ...
    'nonlinear=',20,'maxIter=',5000,'ignoreShocks=',true, ...
    'plan=',p); %?plan?

[s3.ey,s4.ey]

max(abs(s3.ey - s4.ey))
max(abs(s3.epi - s4.epi))
max(abs(s3.er - s4.er))

%% Save Everything for Further Use
%
% Save the model object and the simulated databases to a mat (binary) file
% for further use in other files.

save stochastic_simulations m s1 s2 s3;

%% Help on IRIS Functions Used in This File
%
% Use either `help` to display help in the command window, or `idoc`
% to display help in an HTML browser window.
%
%    help model
%    help model/subsasgn
%    help model/sstatedb
%    help model/simulate
%    help plan
%    help plan/plan
%    help plan/exogenise
%    help plan/endogenise
%    help dbase/dbplot
%    help grfun/bottomlegend
%    help grfun/ftitle
